This lab provides a fairly brief introduction to using R for Bayesian analysis. The first section will cover converting some of the models we have looked at in previous labs to Bayesian models using rstan (an interface to the stan language) and INLA. We’ll then go on to look at using INLA for spatial regression models.
R has a large number of add-on packages for Bayesian analysis listed here. Some notable ones are:
The INLA package is not available through the usual CRAN repositories. Instead, you’ll need to install this manually. The following code will install the current stable version of INLA:
install.packages("INLA",repos=c(getOption("repos"),INLA="https://inla.r-inla-download.org/R/stable"), dep=TRUE)
Once this is installed, we can load the libraries that we’ll need for this lab
library(dplyr)
library(ggplot2)
library(lme4)
library(sf)
library(INLA)
library(spdep)
library(tmap)
gap <- read.csv("../datafiles/gapminderData5.csv")
gap$year2 <- gap$year - 1952
gap$gdpPercap2 <- log(gap$gdpPercap)
fit_lm <- lm(lifeExp ~ gdpPercap2, gap)
summary(fit_lm)
##
## Call:
## lm(formula = lifeExp ~ gdpPercap2, data = gap)
##
## Residuals:
## Min 1Q Median 3Q Max
## -32.778 -4.204 1.212 4.658 19.285
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -9.1009 1.2277 -7.413 1.93e-13 ***
## gdpPercap2 8.4051 0.1488 56.500 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.62 on 1702 degrees of freedom
## Multiple R-squared: 0.6522, Adjusted R-squared: 0.652
## F-statistic: 3192 on 1 and 1702 DF, p-value: < 2.2e-16
inla_lm <- inla(lifeExp ~ gdpPercap2,
data = gap)
inla_lm <- inla(lifeExp ~ gdpPercap2,
data = gap,
control.compute = list(dic = TRUE, waic = TRUE),
control.predictor = list(compute = TRUE))
summary(inla_lm)
##
## Call:
## c("inla(formula = lifeExp ~ gdpPercap2, data = gap, control.compute =
## list(dic = TRUE, ", " waic = TRUE), control.predictor = list(compute =
## TRUE))" )
## Time used:
## Pre = 4.79, Running = 0.549, Post = 0.53, Total = 5.87
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) -9.099 1.228 -11.510 -9.099 -6.691 -9.099 0
## gdpPercap2 8.405 0.149 8.113 8.405 8.697 8.405 0
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.017 0.001 0.016 0.017
## 0.975quant mode
## Precision for the Gaussian observations 0.018 0.017
##
## Expected number of effective parameters(stdev): 2.00(0.00)
## Number of equivalent replicates : 851.93
##
## Deviance Information Criterion (DIC) ...............: 11760.48
## Deviance Information Criterion (DIC, saturated) ....: 1707.78
## Effective number of parameters .....................: 3.03
##
## Watanabe-Akaike information criterion (WAIC) ...: 11760.97
## Effective number of parameters .................: 3.51
##
## Marginal log-Likelihood: -5899.79
## Posterior marginals for the linear predictor and
## the fitted values are computed
In the linear model, we got an estimate of the residual standard deviation of around 7.6. This is a measure of the size of the errors from the model. In the INLA output, we get an estimate of the residual precision instead, which is the inverse of the variance. To compare, we can invert this and take the squre root to convert back to standard deviation units, and to check that the models are estimating the same value:
sqrt(1 / 0.017)
## [1] 7.66965
plot(inla_lm, plot.fixed.effects = TRUE,
plot.random.effects = FALSE,
plot.hyperparameters = FALSE,
plot.predictor = FALSE)
plot(inla_lm, plot.fixed.effects = FALSE,
plot.random.effects = FALSE,
plot.hyperparameters = TRUE,
plot.predictor = FALSE)
plot(gap$lifeExp, inla_lm$summary.linear.predictor$mean)
abline(0,1)
## Binomial regression
irished <- read.csv("../datafiles/irished.csv")
irished$DVRT.cen <- irished$DVRT - mean(irished$DVRT)
fit_glm <- glm(lvcert ~ DVRT.cen,
data = irished,
family = binomial(link = 'logit'))
summary(fit_glm)
##
## Call:
## glm(formula = lvcert ~ DVRT.cen, family = binomial(link = "logit"),
## data = irished)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.9077 -0.9810 -0.4543 1.0307 2.1552
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.278422 0.099665 -2.794 0.00521 **
## DVRT.cen 0.064369 0.007576 8.496 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 686.86 on 499 degrees of freedom
## Residual deviance: 593.77 on 498 degrees of freedom
## AIC: 597.77
##
## Number of Fisher Scoring iterations: 3
inla_glm <- inla(lvcert ~ DVRT.cen,
data = irished,
family = "binomial",
control.compute = list(dic = TRUE, waic = TRUE),
control.predictor = list(compute = TRUE))
summary(inla_glm)
##
## Call:
## c("inla(formula = lvcert ~ DVRT.cen, family = \"binomial\", data =
## irished, ", " control.compute = list(dic = TRUE, waic = TRUE),
## control.predictor = list(compute = TRUE))" )
## Time used:
## Pre = 5.54, Running = 0.323, Post = 0.383, Total = 6.24
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) -0.279 0.100 -0.476 -0.279 -0.084 -0.278 0
## DVRT.cen 0.065 0.008 0.050 0.064 0.080 0.064 0
##
## Expected number of effective parameters(stdev): 2.00(0.00)
## Number of equivalent replicates : 249.92
##
## Deviance Information Criterion (DIC) ...............: 597.75
## Deviance Information Criterion (DIC, saturated) ....: 597.75
## Effective number of parameters .....................: 1.99
##
## Watanabe-Akaike information criterion (WAIC) ...: 597.79
## Effective number of parameters .................: 2.02
##
## Marginal log-Likelihood: -306.61
## Posterior marginals for the linear predictor and
## the fitted values are computed
plot(inla_glm, plot.fixed.effects = TRUE,
plot.random.effects = FALSE,
plot.hyperparameters = FALSE,
plot.predictor = FALSE)
fit_lmer <- lmer(lifeExp ~ year2 + (1 | country), gap)
summary(fit_lmer)
## Linear mixed model fit by REML ['lmerMod']
## Formula: lifeExp ~ year2 + (1 | country)
## Data: gap
##
## REML criterion at convergence: 9866.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.1692 -0.4910 0.1216 0.6030 2.4194
##
## Random effects:
## Groups Name Variance Std.Dev.
## country (Intercept) 123.15 11.097
## Residual 12.85 3.584
## Number of obs: 1704, groups: country, 142
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 50.51208 0.94547 53.43
## year2 0.32590 0.00503 64.79
##
## Correlation of Fixed Effects:
## (Intr)
## year2 -0.146
inla_lmer <- inla(lifeExp ~ year2 + f(country, model = "iid"),
data = gap)
summary(inla_lmer)
##
## Call:
## "inla(formula = lifeExp ~ year2 + f(country, model = \"iid\"), data =
## gap)"
## Time used:
## Pre = 6.04, Running = 0.558, Post = 0.343, Total = 6.94
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 50.512 0.947 48.651 50.512 52.372 50.512 0
## year2 0.326 0.005 0.316 0.326 0.336 0.326 0
##
## Random effects:
## Name Model
## country IID model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.078 0.003 0.072 0.078
## Precision for country 0.008 0.001 0.006 0.008
## 0.975quant mode
## Precision for the Gaussian observations 0.083 0.078
## Precision for country 0.010 0.008
##
## Expected number of effective parameters(stdev): 141.77(0.154)
## Number of equivalent replicates : 12.02
##
## Marginal log-Likelihood: -4968.37
plot(inla_lmer, plot.fixed.effects = TRUE,
plot.random.effects = FALSE,
plot.hyperparameters = FALSE,
plot.predictor = FALSE)
plot(inla_lmer, plot.fixed.effects = FALSE,
plot.random.effects = TRUE,
plot.hyperparameters = FALSE,
plot.predictor = FALSE, cex = 1)
plot(inla_lmer, plot.fixed.effects = FALSE,
plot.random.effects = FALSE,
plot.hyperparameters = TRUE,
plot.predictor = FALSE,
plot.prior = TRUE)
Read in the data
ufos <- st_read("../datafiles/nuforc/ufo_sf_annual.shp")
## Reading layer `ufo_sf_annual' from data source
## `/Users/u0784726/Dropbox/Data/devtools/geog6000/datafiles/nuforc/ufo_sf_annual.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1421 features and 44 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -124.71 ymin: 24.54235 xmax: -66.98702 ymax: 49.36967
## Geodetic CRS: WGS 84
Select only 2014
ufos_2014 <- ufos %>%
filter(year == 2014)
p1 <- tm_shape(ufos_2014) + tm_fill("count", style = "jenks") +
tm_borders() +
tm_layout(main.title = "UFO sightings 2014")
p2 <- tm_shape(ufos_2014) + tm_fill("population", style = "jenks") +
tm_borders() +
tm_layout(main.title = "State population 2014")
tmap_arrange(p1, p2)
ufo_lm <- lm(count ~ population, ufos_2014)
summary(ufo_lm)
##
## Call:
## lm(formula = count ~ population, data = ufos_2014)
##
## Residuals:
## Min 1Q Median 3Q Max
## -261.79 -30.40 -12.68 22.04 219.68
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 34.524562 14.802688 2.332 0.024 *
## population 0.019879 0.001543 12.882 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 76.66 on 47 degrees of freedom
## Multiple R-squared: 0.7793, Adjusted R-squared: 0.7746
## F-statistic: 165.9 on 1 and 47 DF, p-value: < 2.2e-16
ufo_inla_lm <- inla(count ~ population,
data = ufos_2014,
control.compute = list(dic = TRUE, waic = TRUE, cpo = TRUE),
control.predictor = list(compute = TRUE))
summary(ufo_inla_lm)
##
## Call:
## c("inla(formula = count ~ population, data = ufos_2014, control.compute
## = list(dic = TRUE, ", " waic = TRUE, cpo = TRUE), control.predictor =
## list(compute = TRUE))" )
## Time used:
## Pre = 5.03, Running = 0.313, Post = 0.351, Total = 5.7
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 34.525 14.936 5.058 34.524 63.964 34.525 0
## population 0.020 0.002 0.017 0.020 0.023 0.020 0
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.00 0.00 0.00 0.00
## 0.975quant mode
## Precision for the Gaussian observations 0.00 0.00
##
## Expected number of effective parameters(stdev): 2.00(0.00)
## Number of equivalent replicates : 24.50
##
## Deviance Information Criterion (DIC) ...............: 568.69
## Deviance Information Criterion (DIC, saturated) ....: 58.39
## Effective number of parameters .....................: 3.16
##
## Watanabe-Akaike information criterion (WAIC) ...: 574.64
## Effective number of parameters .................: 7.59
##
## Marginal log-Likelihood: -307.59
## CPO and PIT are computed
##
## Posterior marginals for the linear predictor and
## the fitted values are computed
ufo_inla_lm$summary.fixed
## mean sd 0.025quant 0.5quant 0.975quant
## (Intercept) 34.52460098 14.935909271 5.05803232 34.52414369 63.96379416
## population 0.01987854 0.001557171 0.01680645 0.01987849 0.02294778
## mode kld
## (Intercept) 34.52460098 1.497393e-05
## population 0.01987854 1.500443e-05
Posterior distribution
ufo_inla_lm$marginals.fixed
ufo_inla_lm$marginals.hyperpar
Plot posterior distribution of intercept
plot_df <- as.data.frame(inla.tmarginal(function(x) x,
ufo_inla_lm$marginals.fixed$`(Intercept)`))
ggplot(plot_df, aes(x = x, y = y)) +
geom_line(size = 2) +
scale_x_continuous("Intercept") +
scale_y_continuous("P") +
theme_bw()
Plot posterior distribution of population coefficient
plot_df <- as.data.frame(inla.tmarginal(function(x) x,
ufo_inla_lm$marginals.fixed$population))
ggplot(plot_df, aes(x = x, y = y)) +
geom_line(size = 2) +
scale_x_continuous("Population") +
scale_y_continuous("P") +
theme_bw()
Plot posterior distribution of residual variance
plot_df <- as.data.frame(inla.tmarginal(function(x) 1/x,
ufo_inla_lm$marginals.hyperpar$`Precision for the Gaussian observations`))
ggplot(plot_df, aes(x = x, y = y)) +
geom_line(size = 2) +
scale_x_continuous("Residual variance") +
scale_y_continuous("P") +
theme_bw()
Add an informed prior on slope
ufo_inla_lm2 <- inla(count ~ population,
data = ufos_2014,
control.fixed = list(mean = 0, prec = 0.2,
mean.intercept = 0, prec.intercept = 0.2),
control.family=list(hyper=list(prec=list(prior='loggamma',
param=c(0.1,0.1)))),
control.compute = list(dic = TRUE, waic = TRUE),
control.predictor = list(compute = TRUE))
summary(ufo_inla_lm2)
##
## Call:
## c("inla(formula = count ~ population, data = ufos_2014, control.compute
## = list(dic = TRUE, ", " waic = TRUE), control.predictor = list(compute
## = TRUE), control.family = list(hyper = list(prec = list(prior =
## \"loggamma\", ", " param = c(0.1, 0.1)))), control.fixed = list(mean =
## 0, prec = 0.2, ", " mean.intercept = 0, prec.intercept = 0.2))")
## Time used:
## Pre = 5.03, Running = 0.259, Post = 0.346, Total = 5.63
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 0.713 2.217 -3.642 0.713 5.062 0.713 0
## population 0.022 0.001 0.020 0.022 0.025 0.022 0
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.00 0.00 0.00 0.00
## 0.975quant mode
## Precision for the Gaussian observations 0.00 0.00
##
## Expected number of effective parameters(stdev): 1.02(0.004)
## Number of equivalent replicates : 48.01
##
## Deviance Information Criterion (DIC) ...............: 571.62
## Deviance Information Criterion (DIC, saturated) ....: 51.66
## Effective number of parameters .....................: 2.08
##
## Watanabe-Akaike information criterion (WAIC) ...: 577.02
## Effective number of parameters .................: 6.33
##
## Marginal log-Likelihood: -295.48
## Posterior marginals for the linear predictor and
## the fitted values are computed
ufos_2014$idarea <- 1:nrow(ufos_2014)
ufos_nb <- poly2nb(ufos_2014)
ufos_W <- nb2mat(ufos_nb)
plot(st_geometry(ufos_2014), reset = FALSE)
plot(ufos_nb, st_coordinates(st_centroid(ufos_2014)), add = TRUE, col = 2, lwd = 2)
ufo_inla_bym <- inla(count ~ population + f(idarea, model = "bym2", graph = ufos_W),
data = ufos_2014,
control.compute = list(dic = TRUE, waic = TRUE),
control.predictor = list(compute = TRUE)
)
summary(ufo_inla_bym)
##
## Call:
## c("inla(formula = count ~ population + f(idarea, model = \"bym2\", ", "
## graph = ufos_W), data = ufos_2014, control.compute = list(dic = TRUE,
## ", " waic = TRUE), control.predictor = list(compute = TRUE))" )
## Time used:
## Pre = 28.2, Running = 0.478, Post = 0.402, Total = 29.1
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 40.013 6.129 27.891 40.041 51.97 40.097 0
## population 0.019 0.001 0.018 0.019 0.02 0.019 0
##
## Random effects:
## Name Model
## idarea BYM2 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 3.08e+04 1.48e+04 9292.433 2.86e+04
## Precision for idarea 1.00e-03 0.00e+00 0.001 1.00e-03
## Phi for idarea 3.54e-01 1.03e-01 0.188 3.41e-01
## 0.975quant mode
## Precision for the Gaussian observations 6.56e+04 2.28e+04
## Precision for idarea 1.00e-03 1.00e-03
## Phi for idarea 5.84e-01 3.05e-01
##
## Expected number of effective parameters(stdev): 49.00(0.00)
## Number of equivalent replicates : 1.00
##
## Deviance Information Criterion (DIC) ...............: -295.17
## Deviance Information Criterion (DIC, saturated) ....: 125.63
## Effective number of parameters .....................: 62.82
##
## Watanabe-Akaike information criterion (WAIC) ...: -318.00
## Effective number of parameters .................: 30.41
##
## Marginal log-Likelihood: -463.53
## Posterior marginals for the linear predictor and
## the fitted values are computed
Add a prior on the spatial term.
prior_bym2 <- list(
prec = list(
prior = "pc.prec",
param = c(1, 0.01)),
phi = list(
prior = "pc",
param = c(0.5, 2 / 3))
)
ufo_inla_bym2 <- inla(count ~ population + f(idarea, model = "bym2",
graph = ufos_W,
hyper = prior_bym2),
data = ufos_2014,
control.compute = list(dic = TRUE, waic = TRUE),
control.predictor = list(compute = TRUE),
control.inla = list(strategy = "adaptive",int.strategy = "eb")
)
summary(ufo_inla_bym2)
##
## Call:
## c("inla(formula = count ~ population + f(idarea, model = \"bym2\", ", "
## graph = ufos_W, hyper = prior_bym2), data = ufos_2014, control.compute
## = list(dic = TRUE, ", " waic = TRUE), control.predictor = list(compute
## = TRUE), control.inla = list(strategy = \"adaptive\", ", " int.strategy
## = \"eb\"))")
## Time used:
## Pre = 25.8, Running = 0.461, Post = 0.368, Total = 26.6
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 34.525 14.124 6.794 34.524 62.233 34.525 0
## population 0.020 0.001 0.017 0.020 0.023 0.020 0
##
## Random effects:
## Name Model
## idarea BYM2 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.000 0.000 0.000 0.000
## Precision for idarea 1.192 0.413 0.587 1.124
## Phi for idarea 0.048 0.016 0.022 0.047
## 0.975quant mode
## Precision for the Gaussian observations 0.000 0.000
## Precision for idarea 2.187 1.001
## Phi for idarea 0.082 0.044
##
## Expected number of effective parameters(stdev): 2.01(0.00)
## Number of equivalent replicates : 24.40
##
## Deviance Information Criterion (DIC) ...............: 566.34
## Deviance Information Criterion (DIC, saturated) ....: 55.66
## Effective number of parameters .....................: 2.01
##
## Watanabe-Akaike information criterion (WAIC) ...: 569.79
## Effective number of parameters .................: 4.78
##
## Marginal log-Likelihood: -283.33
## Posterior marginals for the linear predictor and
## the fitted values are computed
Plot spatial effect
ufos_2014$u <- ufo_inla_bym2$summary.random$idarea$mean[1:49]
tm_shape(ufos_2014) + tm_fill("u", style = "jenks") +
tm_borders() +
tm_layout(main.title = "State population 2014")
Plot Bayesian estimates of sightings
ufos_2014$yhat <- ufo_inla_bym2$summary.linear.predictor$`0.025quant`
tm_shape(ufos_2014) + tm_fill("yhat", style = "jenks") +
tm_borders() +
tm_layout(main.title = "State population 2014")
ufos_2014$Ei <- (sum(ufos_2014$count) / sum(ufos_2014$population)) * ufos_2014$population
ufos_2014$SIR <- ufos_2014$count / ufos_2014$Ei
tm_shape(ufos_2014) + tm_fill("SIR") +
tm_borders() + tm_layout(main.title = "UFO SIR 2014")
ufo_inla_rr <- inla(count ~ 1 + f(idarea, model = "bym2",
graph = ufos_W,
hyper = prior_bym2),
data = ufos_2014, family = "poisson", E = Ei,
control.compute = list(dic = TRUE, waic = TRUE),
control.predictor = list(compute = TRUE),
control.inla = list(strategy = "adaptive",int.strategy = "eb")
)
summary(ufo_inla_rr)
##
## Call:
## c("inla(formula = count ~ 1 + f(idarea, model = \"bym2\", graph =
## ufos_W, ", " hyper = prior_bym2), family = \"poisson\", data =
## ufos_2014, ", " E = Ei, control.compute = list(dic = TRUE, waic =
## TRUE), ", " control.predictor = list(compute = TRUE), control.inla =
## list(strategy = \"adaptive\", ", " int.strategy = \"eb\"))")
## Time used:
## Pre = 25.6, Running = 0.379, Post = 0.36, Total = 26.3
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 0.033 0.047 -0.06 0.033 0.126 0.033 0
##
## Random effects:
## Name Model
## idarea BYM2 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Precision for idarea 4.759 1.30 2.660 4.610 7.738 4.328
## Phi for idarea 0.568 0.24 0.115 0.584 0.949 0.771
##
## Expected number of effective parameters(stdev): 45.04(0.00)
## Number of equivalent replicates : 1.09
##
## Deviance Information Criterion (DIC) ...............: 420.40
## Deviance Information Criterion (DIC, saturated) ....: 106.46
## Effective number of parameters .....................: 45.18
##
## Watanabe-Akaike information criterion (WAIC) ...: 414.14
## Effective number of parameters .................: 28.14
##
## Marginal log-Likelihood: -236.32
## Posterior marginals for the linear predictor and
## the fitted values are computed
ufos_2014$RR <- ufo_inla_rr$summary.fitted.values[, "mean"]
tm_shape(ufos_2014) + tm_fill("RR", palette = "-inferno") +
tm_borders() +
tm_layout(main.title = "RR")
ufos_2014$exc <- sapply(ufo_inla_rr$marginals.fitted.values,
function(marg){1 - inla.pmarginal(q = 2, marginal = marg)})
tm_shape(ufos_2014) + tm_fill("exc", palette = "Blues") +
tm_borders() +
tm_layout(main.title = "Excedence probability (RR > 2)")
airports <- read.csv("../datafiles/nuforc/airports.csv")
airports$lairports <- log(airports$airports)
ufos_2014 <- merge( ufos_2014, airports, by.x = "postal", by.y = "state")
ggplot(ufos_2014, aes(x = airports, y = count)) +
scale_x_log10("Airports") + scale_y_log10("Sightings") +
geom_point()
prior_bym2 <- list(
prec = list(
prior = "pc.prec",
param = c(0.6 / 0.31, 0.01)),
phi = list(
prior = "pc",
param = c(0.005, 1 / 10))
)
ufo_inla_rr2 <- inla(count ~ lairports +
f(idarea, model = "bym2",
graph = ufos_W,
hyper = prior_bym2),
data = ufos_2014, family = "poisson", E = Ei,
control.compute = list(dic = TRUE, waic = TRUE),
control.predictor = list(compute = TRUE),
control.inla = list(strategy = "adaptive",int.strategy = "eb")
)
summary(ufo_inla_rr2)
##
## Call:
## c("inla(formula = count ~ lairports + f(idarea, model = \"bym2\", ", "
## graph = ufos_W, hyper = prior_bym2), family = \"poisson\", ", " data =
## ufos_2014, E = Ei, control.compute = list(dic = TRUE, ", " waic =
## TRUE), control.predictor = list(compute = TRUE), ", " control.inla =
## list(strategy = \"adaptive\", int.strategy = \"eb\"))" )
## Time used:
## Pre = 23.2, Running = 0.344, Post = 0.362, Total = 23.9
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 0.419 0.482 -0.529 0.420 1.364 0.421 0
## lairports -0.068 0.084 -0.233 -0.068 0.097 -0.068 0
##
## Random effects:
## Name Model
## idarea BYM2 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Precision for idarea 4.915 1.212 2.915 4.792 7.638 4.56
## Phi for idarea 0.142 0.144 0.005 0.093 0.542 0.01
##
## Expected number of effective parameters(stdev): 45.72(0.00)
## Number of equivalent replicates : 1.07
##
## Deviance Information Criterion (DIC) ...............: 420.62
## Deviance Information Criterion (DIC, saturated) ....: 106.67
## Effective number of parameters .....................: 45.87
##
## Watanabe-Akaike information criterion (WAIC) ...: 414.30
## Effective number of parameters .................: 28.51
##
## Marginal log-Likelihood: -242.86
## Posterior marginals for the linear predictor and
## the fitted values are computed